//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   BuiltInCrossSections.cc
 * \author Jeremy Roberts
 * \date   Jul 17, 2011
 * \brief  
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <vector>
#include <iostream>

#include "BuiltInCrossSections.hh"
#include "../utilities/DBC.hh"
#include "../utilities/SP.hh"

namespace slabtran
{

// Constructor
BuiltInCrossSections::BuiltInCrossSections(const int setid)
{
  Require(setid >= 0 && setid < END_BUILTINXS);
  // point to the right built in data
  const int *data;
  const double *sigma;
  const string *name;
  switch (setid)
  {
  case XS1G1M:
    data = xs1g1m;
    sigma = xs1g1m_sigma;
    name = xs1g1m_name;
    d_mat_name = "xs1g1m";
    break;
  case XS1G2M:
    data = xs1g2m;
    sigma = xs1g2m_sigma;
    name = xs1g2m_name;
    d_mat_name = "xs1g2m";
    break;
  case XS2G4M:
    data = xs2g4m;
    sigma = xs2g4m_sigma;
    name = xs2g4m_name;
    d_mat_name = "xs2g4m";
    break;
  case XS6G1M:
    data = xs6g1m;
    sigma = xs6g1m_sigma;
    name = xs6g1m_name;
    d_mat_name = "xs6g1m";
    break;
  default:
    throw util::assertion("Invalid built in cross-section id.");
  }
  set_data(data,sigma,name);
}

//
BuiltInCrossSections::SP_materials BuiltInCrossSections::build()
{
  // build the actual library
  SP_materials mat;
  mat = new Materials(d_number_materials,d_mat_name);
  Materials::SP_XS y; // place holder for element of mat

  //-------------------------------------------------------------------------
  // build the cross-sections for each material
  for (int m = 0; m < d_number_materials; m++)
  {
    //-------------------------------------------------------------------------
    // create a new cross-section object
    //mat->get_p(m) = new Cross_Section(d_number_groups, d_legendre_order);
    //Cross_Section y = mat->get_p(m);
    y = new Cross_Section(d_number_groups, d_legendre_order);
    //y = mat->get_p(m); // assign placeholder

    //-------------------------------------------------------------------------
    // Name it
    y->set_name(d_name[m]);

    //-------------------------------------------------------------------------
    // fill group data
    int col; // used below as a column placeholder
    for (int g = 0; g < d_number_groups; g++)
    {
      col = 0;
      int row = g + m * d_number_groups;
      // assign group data
      y->sigma(g)    = d_sigma[row][col++];
      y->sigma_a(g)  = d_sigma[row][col++];
      y->sigma_f(g)  = d_sigma[row][col++];
      y->nu(g)       = d_sigma[row][col++];
      y->chi(g)      = d_sigma[row][col++];
    }

    //-------------------------------------------------------------------------
    // assign lower and upper bounds by:
    //  -- loop over gprime=0..ng for a g to get first nonzero = d_group_lower_bound[0][g]
    //  -- loop over gprime=ng..0 ...                          = d_group_upper_bound[0][g]
    //  -- loop over g=0..ng
    //  -- loop over g=ng..0
    //
    //                               s0<0  s0<1  s0<2  s0<3  upscatter
    // e.g. for one material:        s1<0  s1<1  s1<2  s1<3
    //                               s2<0  s2<1  s2<2  s2<3
    //                  downscatter  s3<0  s3<1  s3<2  s3<3
    // supp
    for (int g = 0; g < d_number_groups; g++)                  // given row
    {
      int row = g + m * d_number_groups;
      for (int gp = 0; gp < d_number_groups; gp++)             // going right
      {
        // must have a non-negative zeroth moment
        Check(d_sigma[row][col+gp*(d_legendre_order+1)] >= 0.0);
        if (d_sigma[row][col+gp*(d_legendre_order+1)] > 0.0)
        {
          y->group_lower_bound(g, 0) = gp;                      // keep lowest nonzero group (that downscatters to me)
          break;
        }
      }
      for (int gp = d_number_groups-1; gp >= 0; gp--)          // going left
      {
        // must have a non-negative zeroth moment
        Check(d_sigma[row][col+gp*(d_legendre_order+1)] >= 0.0);
        if (d_sigma[row][col+gp*(d_legendre_order+1)] > 0.0)
        {
          y->group_upper_bound(g, 0) = gp;                      // keep highest nonzero group (that upscatters to me)
          break;
        }
      }
    }
    for (int gp = 0; gp < d_number_groups; gp++)                // given column
    {
      for (int g = 0; g < d_number_groups; g++)                // going down
      {
        int row = g+m*d_number_groups;
        if (d_sigma[row][col+gp*(d_legendre_order+1)] > 0.0)
        { // note the second argument to group_upper_bound
          y->group_lower_bound(gp, 1) = g;                      // keep lowest nonzero group (to which i upscatter)
          break;
        }
      }
      for (int g = d_number_groups-1; g >= 0; g--)             // going up
      {
        int row = g+m*d_number_groups;
        if (d_sigma[row][col+gp*(d_legendre_order+1)] > 0.0)
        { // note the second argument to group_upper_bound
          y->group_upper_bound(gp, 1) = g;                      // keep highest nonzero group (to which i downscatter)
          break;
        }
      }
    }

    //-------------------------------------------------------------------------
    // build S rows and place in nonzero elements
    vector<double> block;
    for (int g = 0; g < d_number_groups; g++) // all rows
    {
      int row = g+m*d_number_groups;
      for (int gp = y->group_lower_bound(g,0); gp <= y->group_upper_bound(g,0); gp++)
      {
        block.clear();                        // clear and
        block.resize(d_legendre_order + 1); // allocate
        // row for g has g<-1, g<-2, ... g<-G
        for (int l = 0; l < d_legendre_order + 1; l++)
        {
          block[l] = d_sigma[row][col+gp*(d_legendre_order+1)+l];
        }
        // insert block
        y->sigma_s(g,gp) = block;
      }
    }
    //delete y;
    // point mat element to y
    mat->get_p(m) = y;

  } // end material loop

  // return the material library
  return mat;
} // end build

// PRIVATE -- IMPLEMENTATION

// fill data from a constant set
void BuiltInCrossSections::set_data(const int *data, const double *sigma, const string *name)
{
  // set integer data
  d_number_groups = data[0];        //std::cout << " d_number_groups = " << d_number_groups << std::endl;
  d_number_materials = data[1];     //std::cout << " d_number_materials = " << d_number_materials << std::endl;
  d_legendre_order = data[2];       //std::cout << " d_legendre_order = " << d_legendre_order << std::endl;
  // resize name and fill
  d_name.resize(d_number_materials);
  for (int m = 0; m < d_number_materials; m++)
    d_name[m] = name[m];
  // resize sigma
  d_sigma.resize(d_number_groups*d_number_materials);
  int col = 5+d_number_groups*(d_legendre_order+1);
  for (int i = 0; i < d_sigma.size(); i++)
  {
    // T,A,F,nu,chi + scattering
    d_sigma[i].resize(col);
    for (int j = 0; j < d_sigma[0].size(); j++)
    {
      d_sigma[i][j] = sigma[i*col + j]; //std::cout  << d_sigma[i][j] << " ";
    }
    //std::cout << std::endl;
  }
}

} // end namespace slabtran


//---------------------------------------------------------------------------//
//              end of BuiltInCrossSections.cc
//---------------------------------------------------------------------------//
